home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / polywarp.pro < prev    next >
Encoding:
Text File  |  1997-07-08  |  3.3 KB  |  134 lines
  1. ; $Id: polywarp.pro,v 1.3 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1983-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5.  
  6. pro POLYWARP, XI, YI, XO, YO, DEGREE, KX, KY
  7. ;+
  8. ; NAME:
  9. ;    POLYWARP
  10. ;
  11. ; PURPOSE:
  12. ;    Perform polynomial spatial warping.
  13. ;
  14. ;    Using least squares estimation, determine the coefficients Kx[i,j] 
  15. ;    and Ky[i,j] of the polynomial functions:
  16. ;        Xi = sum over i and j of:  Kx[i,j] * Xo^j * Yo^i
  17. ;        Yi = sum over i and j of:  Ky[i,j] * Xo^j * Yo^i
  18. ;
  19. ;    Kx and Ky can be used as inputs P and Q to the built-in function
  20. ;    POLY_2D.
  21. ;    
  22. ; CATEGORY:
  23. ;    Image processing.
  24. ;
  25. ; CALLING SEQUENCE:
  26. ;    POLYWARP, Xi, Yi, Xo, Yo, Degree, Kx, Ky
  27. ;
  28. ; INPUTS:
  29. ;    Xi, Yi:    The vectors of x,y coordinates to be fit as a function 
  30. ;        of Xo and Yo.
  31. ;
  32. ;    Xo, Yo:    The vectors of x,y independent coordinates.  These vectors 
  33. ;        must have the same number of elements as Xi and Yi.
  34. ;
  35. ;    Degree:    The degree of the fit.  The number of coordinate pairs must be
  36. ;        greater than or equal to (Degree+1)^2.
  37. ;
  38. ; OUTPUTS:
  39. ;    Kx:    The array of coefficients for Xi as a function of (xo,yo).
  40. ;        This parameter is returned as a (Degree+1) by (Degree+1)
  41. ;        element array.
  42. ;
  43. ;    Ky:    The array of coefficients for yi.  This parameter is returned
  44. ;        as a (Degree+1) by (Degree+1) element array.
  45. ;
  46. ; COMMON BLOCKS:
  47. ;    None.
  48. ;
  49. ; SIDE EFFECTS:
  50. ;    None.
  51. ;
  52. ; RESTRICTIONS:
  53. ;    None.
  54. ;
  55. ; PROCEDURE:
  56. ;    See:    Computer Image Processing and Recognition, Ernest L. Hall,
  57. ;        Academic Press, 1979, Pages 186-188.
  58. ;
  59. ;    Xi and Yi are expressed as polynomials of Xo, Yo:
  60. ;        Xi = Kx[i,j] * Xo^j * Yo^i   Summed for i,j = 0 to degree.
  61. ;    And
  62. ;        Yi = Ky[i,j] * Xo^j * Yo^i.
  63. ;
  64. ;    This coordinate transformation may be then used to
  65. ;    map from Xo, Yo coordinates into Xi, Yi coordinates.
  66. ;
  67. ; EXAMPLE:
  68. ;    The following example shows how to display an image and warp it
  69. ;    using the POLYWARP and POLY_2D routines.
  70. ;
  71. ;    Create and display the original image by entering:
  72. ;
  73. ;        A = BYTSCL(SIN(DIST(250)))
  74. ;        TVSCL, A
  75. ;
  76. ;    Now set up the Xi's and Yi's.  Enter:
  77. ;
  78. ;        XI = [24, 35, 102, 92]
  79. ;        YI = [81, 24, 25, 92]
  80. ;
  81. ;    Enter the Xo's and Yo's:
  82. ;
  83. ;        XO = [61, 62, 143, 133]
  84. ;        YO = [89, 34, 38, 105]
  85. ;
  86. ;    Run POLYWARP to obtain a Kx and Ky:
  87. ;
  88. ;        POLYWARP, XI, YI, XO, YO, 1, KX, KY
  89. ;
  90. ;    Create a warped image based on Kx and Ky with POLY_2D:
  91. ;
  92. ;        B = POLY_2D(A, KX, KY)
  93. ;
  94. ;    Display the new image:
  95. ;
  96. ;        TV, B
  97. ;
  98. ; MODIFICATION HISTORY:
  99. ;    DMS, RSI, Dec, 1983.
  100. ;-
  101. ;
  102. on_error,2                      ;Return to caller if an error occurs
  103. m = n_elements(xi)        ;# of points..
  104. if (m ne n_elements(yi)) or (n_elements(xo) ne n_elements(yo)) $
  105.     or (m ne n_elements(xo)) then begin
  106.         message,'Inconsistent number of elements.'
  107.         endif
  108. ;
  109. n = degree        ;use halls notation
  110. n2=(n+1)^2
  111. if n2 gt m then message, '# of points must be ge (degree+1)^2.'
  112. ;
  113. x = dblarr(2,m)        ;x array
  114. u = x
  115. x = double([transpose(xi[*]),transpose(yi[*])])
  116. u = double([transpose(xo[*]),transpose(yo[*])])
  117. ;
  118. ut = dblarr(n2,m)    ;transpose of U
  119. u2i = dblarr(n+1)    ;[1,u2i,u2i^2,...]
  120. for i=0L,m-1 do begin
  121.     u2i[0]=1.    ;init u2i
  122.     zz = u[1,i]
  123.     for j=1,n do u2i[j]=u2i[j-1]*zz
  124.     ut[0,i]= u2i    ;evaluate 0 th power separately
  125.     for j=1,n do ut[j*(n+1),i]=u2i*u[0,i]^j ;fill ut=u0i^j * U2i
  126.     endfor
  127. ;
  128. uu = transpose(ut)    ;big u
  129. kk = invert(ut#uu)#ut    ;solve equation
  130. kx = fltarr(n+1,n+1) + float(kk # transpose(x[0,*]))    ;g1, make 2d square
  131. ky = fltarr(n+1,n+1) + float(kk # transpose(x[1,*]))    ;g2
  132. return
  133. end
  134.